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We propose an exactly solvable multisite interaction spin-1/2 Ising-Heisenberg model on a tri- 
angnlated Husimi lattice for the rigorous studies of chaotic entanglement. By making use of the 
generalized star-triangle transformation, we map the initial model onto an effective Ising one on 
a Husimi lattice, which we solve then exactly by applying the recursive method. Expressing the 
entanglement of the Heisenberg spins, that we quantify by means of the concurrence, in terms of 
the magnetic quantities of the system, we demonstrate its bifurcation and chaotic behavior. Fur¬ 
thermore, we show that the underlying chaos may slightly enhance the amount of the entanglement, 
and present on the phase diagram the transition lines from the uniform to periodic and from the 
periodic to chaotic regimes. 
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I. INTRODUCTION 

Quantum entanglement, that stands for the intrinsic 
nonlocal correlations inherent in the quantum theory, 
plays a crucial role in the understanding of the funda¬ 
mentals of quantum physics as well as serves as the key 
resource for implementing practical applications within 
the held of quantum technologies. Specihcally, entan¬ 
glement is vital to diverse secure quantum communica¬ 
tion and teleportation schemes [1,2] and is at the heart 
of quantum computation, quantum metrology, quantum 
imaging [3-5], etc. 

Various protocols, applicable to continuous- and 
discrete-variable systems have been proposed recently, 
allowing one to generate entanglement shared between 
quantum memories [6] and encoded in states of light helds 
[7] or, alternatively, to create systems of strongly entan¬ 
gled superconducting or cavity QED qubits, quantum 
dots, highly excited Rydberg atoms [8-10], etc. Mean¬ 
while, solid-state systems and, particularly, magnetic ma¬ 
terials, are of signihcant importance in this respect as 
they appear to be a source of entanglement too, even 
on a macroscopic level [11, 12]. Furthermore, entangle¬ 
ment turns out to be a powerful tool here for the char¬ 
acterization of quantum phase transitions as it encodes 
all the information shared between subparts of a system 
that may fail to be captured by means of ordinary cor¬ 
relation functions (this can be the case in some exotic 
quantum ground states [13]). In this respect, spin mod¬ 
els are widely used to describe the properties of such 
solid-state systems where intercoupled (e.g., by means of 
the Heisenberg exchange interaction) spins are nested in 
the sites of a specific lattice structure. The entanglement 
features of various spin-lattice models of this kind have 
been thoroughly studied recently [14-17], revealing, par¬ 
ticularly, a strong connection of magnetic properties and 
entanglement [18-20] that allows one to witness the lat¬ 
ter experimentally (e.g., by means of heat capacity and 


magnetic susceptibility measurements [11, 21]). 

Meanwhile, another intriguing question that we ad¬ 
dress in this paper, is the connection of quantum entan¬ 
glement and chaos. As is known, the latter also plays an 
important role in various areas of research, and particu¬ 
larly in the field of quantum information. For instance, 
the uncontrolled interactions between the qubits of a 
quantum computer, being above some critical strength, 
induce chaotic behavior that may break down the effi¬ 
ciency of the computer or bring about exponential decay 
of the entanglement [22, 23]. Consequently, the complex¬ 
ity of error correction codes increases and the correction 
time, however, is independent of whether the system is 
chaotic or not [24]. On the other hand, quantum chaos 
may also ensue in specific quantum protocols, such as 
the Grover’s search algorithm, or quantum Fourier trans¬ 
form, giving rise to a peculiar combination of quantum 
signatures of chaos and integrability [25]. Meanwhile, 
another suitable platform for the combined studies of en¬ 
tanglement and chaotic behavior is nonlinear dissipative 
oscillators [26] and so-called coupled tops [27], where the 
underlying classical chaos may enhance and determine 
the amount of the entanglement shared between the cou¬ 
pled systems [28, 29]. Moreover, quantum signatures of 
chaotic behavior in a kicked top have been demonstrated 
for trapped caesium atoms, showing that if prepared in 
the chaotic sea, the electron and nuclear spins of the 
atoms get rapidly entangled [30]. Finally, some other 
studies, e.g., for iV-atom Jaynes-Cummings model have 
also shown that the entanglement rate is considerably 
enhanced for chaotic initial conditions [31]. 

On our part we propose a specific spin-lattice model to 
study the relation among its chaotic behavior, entangle¬ 
ment, and magnetic properties. Specifically, the struc¬ 
ture that we consider is of recursive nature and consists 
of two non-equivalent types of sites: the nodal ones that 
constitute a Husimi tree, and trimeric units that form 
embedded triangles [see Fig. 1(a)]. Being constructed 
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in such a way, it is thus natural to call this structure 
a triangulated Husimi tree and its deep interior a tri¬ 
angulated Husimi lattice (THL) (we note that a simi¬ 
lar system has been proposed very recently in Ref. [32]). 
Next, taking into account that the rigorous treatment 
of purely Heisenberg models are mostly unattainable, 
we adopt here a spin- 1/2 Ising-Heisenberg model on a 
THL such that the Ising spins are nested on the nodal 
sites of the lattice, while the Heisenberg variables are 
situated on the vertices of the embedded triangles. A 
distinguishable feature of this model is that it can be 
solved exactly by combining the generalized star-triangle 
transformation [33, 34] and the recursive method [35- 
41]. Even more, if three-site interactions between nodal 
Ising spins are included, the model exhibits chaotic be¬ 
havior (this is rather different from the case studied in 
Ref. [32] where only two-site interactions are considered 
in the absence of the external magnetic field). Although 
this fact has been well known for an Ising model on a 
simple Husimi lattice [37, 39, 40], it brings about dras¬ 
tic changes in the properties of the Heisenberg trimer of 
the THL. Namely, at sufficiently low temperatures, its 
magnetization, and even more interestingly, the entan¬ 
glement (that we quantify by means of the concurrence 
[42]) exhibit period-doubling bifurcations, chaotic behav¬ 
ior, and periodic windows. It is worth noting, however, 
that the ground state of the Heisenberg triangle remains 
rigid and that it is the Ising sublattice (to which the 
triangles are coupled) that gives rise to chaos. Further¬ 
more, as is shown below, the spin-spin entanglement of 
the triangle may be slightly enhanced by means of the 
underlying chaotic behavior when compared to the case 
of only three qubits in an external magnetic field. Nev¬ 
ertheless, we point out that quantum correlations cannot 
be induced solely by chaos and that these two are rather 
interconnected: chaos ensues when the entangled W state 
becomes a ground state of the Heisenberg trimer. 

Finally, we note that the Husimi structure, which we 
study here, is a good approximation to more realistic 
models with multisite interactions on ordinary lattices, 
being more reliable than the standard Bethe and mean- 
field approaches [36, 43]. Additionally, it can be used 
for the description of various polymers, and particularly, 
of RNA-like molecules [44] . Meanwhile, triangulated lat¬ 
tices, in their turn, are also of current great interest since 
they may describe real materials (e.g., copper based co¬ 
ordination compounds Cu 9 X 2 (cpa )6 • nH 2 0; X = F, Cl, 
Br [20, 45]), and additionally turn out to be an effi¬ 
cient ground for the understanding of specific mecha¬ 
nisms of quantum ordering at low temperatures [33]. 
Consequently, along with the aforementioned features, 
these arguments reinforce our interest in the studies of 
exactly solvable Ising-Heisenberg models on triangulated 
Husimi structures. 

The paper is organized as follows: in Sec. H we intro¬ 
duce the multisite interaction spin-I/2 Ising-Heisenberg 
model on a THL and show its exact solution by means 
of the generalized star-triangle transformation and the 


recursive approach. In Sec. HI we discuss the magnetic 
and entanglement properties of our model, and reveal 
their chaotic nature. Finally, we draw our conclusions in 
Sec. IV. 


II. MODEL 


We consider a recursive structure that is formed by a 
basic triangular unit bearing Ising variables on its ver¬ 
tices, and Heisenberg-type spins situated on its bonds, 
which, in their turn, form an embedded triangle. Taking 
this building block as the zeroth shell, we construct the 
next one, by connecting q — 1 such units on each site of 
the bigger triangle, and continue this process to develop 
higher-generation shells. Eventually we end up with a 
structure that can be called a triangulated Husimi tree, 
since the Ising spins of this system constitute a regular 
Husimi tree with a coordination number q [Fig. 1(a)]. 
Note that the number of sites both at the surface and 
inside the tree grows exponentially, and for that reason 
we consider the properties of the system in its deep inte¬ 
rior (for more details on this see, e.g.. Ref. [46]), i.e., on 
a triangulated Husimi lattice. 

The definition of the spin-I/2 Ising-Heisenberg model 
on the above structure is now straightforward. The cor¬ 
responding Hamiltonian reads: 
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Here Hk denotes the Hamiltonian of the fcth trian¬ 
gle (IVa being the number of these triangles). Ski = 
{S'^. ,5'^. stands for Heisenberg spin-I/2 operators, 
while for that of the Ising ones (z = I, 2, 3). The first 
term in TLk corresponds to the three-site interaction of a 
strength J 3 between nodal Ising spins, Jm is the strength 
of the Ising-to-Heisenberg coupling, whereas and 
are the interaction strengths of the Heisenberg exchange 
coupling on the {x, y) plane and the z direction, respec¬ 
tively (equivalently, this interaction can be defined as 
of a strength J = with the anisotropy parameter 
A = J^/Jh)- The system is also subject to external 
magnetic fields Hi (note that each Ising spin belongs si¬ 
multaneously to q triangles) and Hu, acting upon the 
Ising- and Heisenberg-type variables, and directed along 
the z axis [see also Fig. 1(a)]. 
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FIG. 1. (Color online) (a) The Ising-Heisenberg model on a triangulated Husimi lattice with coordination number g = 3. The 
red (gray) circles denote the position of the Heisenberg Sk^ = {S^., S ^,, 5^.} spins, and the black ones denote that of the Ising 

/ifc. variables. The full red lines label the Heisenberg-type exchange interaction between Ski spins of the kih triangle, while the 
dashed blue ones stand for the Ising-type coupling of S^i and (i,j = 1,2,3); Sk^fJ-k^ Sk^Hk^ interactions are 

shown only for the central triangle. There is also a three-site exchange coupling that is not shown in this figure, 

(b) The effective Ising model on a Husimi lattice which the above Ising-Heisenberg model is mapped onto. The dashed blue 
lines label the two-site interactions, while the full black arrowed ellipse stands for the three-site interaction (shown only for 
the central triangle). The green dot-dashed lines show how the system is cut apart for employing the recursive method (see 
Sec. IIB for details). 


A. The generalized star-triangle transformation 
and the effective Hamiltonian 


To solve the above defined model exactly we employ 
the generalized star-triangle {Y — A) transformation [32- 
34], that allows one to map the initial Ising-Heisenberg 
model on a THL onto an effective Ising model on a 
Husimi lattice. For that, owing to the fact that the 
triangular Hamiltonians TLk are commutative with one 
another, i.e., ['Hk,'Hi] = 0, we partially factorize the par¬ 
tition function of the initial system to partition functions 
of triangular Hamiltonians (we work in units such that 
the Boltzmann constant is unity): 


*Va-1 

2 =Y.I[ Trfcexp(-Hfc/T) = 
{/**;} ^-0 
Wa-I 

{Mfcj} ^=0 


( 2 ) 


Here T is the temperature, Tr* denotes the trace over 
the degrees of freedom of the three Heisenberg spins of 
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the fcth triangle, whereas the subsequent summation runs 
over all the possible configurations of the nodal Ising 
spins, having the form 


Zk {fJ'ki ) 1**2 ’ “ 


( 3 ) 
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In this equation = ±3/2, ±1/2 are the eigenvalues 
of the total Ising spin operator of the fcth triangle = 
/i|i ^ specific spin configuration, and W (^fc) 


is given as: 
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Having brought the partition function to the form of 
Eq. (3), we apply the generalized star-triangle transfor¬ 
mation, for mapping the initial Ising-Heisenberg model 
onto an effective Ising one: 


= exp 

A exp 
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As the above expression should hold true for any config¬ 
uration of Ising spins, it results in a set of four equations 
(the initial eight ones corresponding to 2^ = 8 configura¬ 
tions of the three nodal spins are reduced to just four due 
to a threefold degeneracy of = 1/2 and = —1/2 
states). Eventually, these equations determine unam¬ 
biguously the mapping parameters through the following 
expressions: 


A = [1E(3/2)1E(-3/2)]1/®[1E(1/2)1F(-1/2)]3/8, 


reff 


Js +rin 
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Consequently, owing to Eq. (3) we rewrite the partition 
function of the initial Hamiltonian of Eq. (1) as 


Z( J3, JiH, , i?H, T, g) = 


which maps the spin-1/2 Ising-Heisenberg model on a 
THE onto a spin-1/2 Ising model on a simple Husimi 
lattice with effective three- and two-site interactions in an 
effective magnetic field. The corresponding Hamiltonian 
reads [see also Fig. 1(b)]: 




A^a-I 

E 


J eff z z z I reff / z z > z z > z 2:\i 
3 MfeiMfcaMfca + J2 (MfeiMfca + + MfeiMfca) + 


^eff 

(^fci + ^J■k2 + Mfea ) 


( 8 ) 


The distinguishable feature of this final model is the 
fact that it can be solved exactly through the recursive 
method, that we present in the next subsection. We also 
make a note here that for the values of parameters that 
we consider below, the effective coupling J|® always re¬ 
mains ferromagnetic. 

Meanwhile, we express all the physical quantities of in¬ 
terest in terms of the above defined effective parameters. 
Specifically, denoting by /a and the free energy per 
triangle of the initial and effective models, from Eq. (7) 
we find: 

fA=ff-T\nA. (9) 


Consequently, the single-site magnetization mi of the 
nodal Ising spins reads (there are 3/g of them in each 
triangle): 


q dfA q 9ff 

3dHi 3 9i7®ff dHi 


( 10 ) 


In other words, mi is equal to the per-cite magnetization 
m®® of the effective model, which is, however, not the case 
for the single-site magnetization mu of the Heisenberg 
sublattice: 


1 dfA 


1 ( Off dff dJf dff dJf\ 

3 dHn ^ dJf OHh ^ dJf dHn J 

T 5 In A _ m®® 5iJ®® e|® dJf ef dJf 

3 dHu q dHu 3 dHu 3 dHu 
T dA 


-f 


( 11 ) 


where ef = and ef = + 

11 ^® three- and two-site spin-spin correlation 
functions of the effective Ising model. 

In what follows we are also interested in the nearest- 
neighbor two-site correlation functions of the Heisenberg 
spins. Note that as the magnetic field is directed along 
the z axis, the correlations in the x and y directions are 
equal. Therefore, the per-pair spin-spin correlation func¬ 
tions read: 


_ 1 5/a _ 

-Syy- ^ - 

1 ( dff an®® dff dJf dff aj|®\ 
3 1^577®® dJ^y ^ dJfdJ^y ^ djf dJ^^y J 
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T dA 
3A dJxy' 


Similarly, the two-site spin-spin correlation function in 
the 2 direction is given as: 
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(13) 
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B. Exact solution of the effective model by means 
of the recursive method 


As mentioned previously, the effective Ising model on 
a Husimi lattice, defined by the Hamiltonian ( 8 ), can be 
solved exactly through the recursive method. The latter 
takes into account the specific — recursive — structure of 
the lattice, that can be cut apart at the central unit (the 
central triangle for the case of a Husimi lattice), being di¬ 
vided into q identical branches [Fig 1(b)]. Consequently, 
the partition function is expressed in terms of the parti¬ 
tion functions of each of those branches [35, 37-40]: 




exp 

MOi ,1^02 'MO 3 


'Jf 

^rM0iM02M03 + 


eff 


^2 
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5r^(Moi)ffr^(M02)5r^(Ai03), 

I 


where the summation runs over the eigenvalues /xg- = 
±1/2 of the Ising spin operators (i.e., over the pos¬ 
sible spin configurations), 5 n(M 0 i) standing for the con¬ 
tribution to the partition function of each of the q — 1 
branches that contain n shells. Next, we repeat the 
cutting-apart procedure once again, and e.g., for g„(/xoi) 
we find: 


ffn(M0i)= X! 
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-j-i. 


Eventually, this procedure results in a recursive relation 
for Xn = 5n(l/2)/g„(-l/2) (see, e.g. Refs. [39, 40, 47]): 
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Although Xn does not have a direct physical meaning, 
the mapping f(x) describes the properties of the sys¬ 
tem. Specifically, the above defined single-site magneti¬ 
zation mi, as well as the three- and two-site correlation 


functions and £ 3 ® of spins situated deep inside the 
tree (where the surface effects are negligible), can be ex¬ 
pressed straightforwardly in terms of x„: 
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In other words, in the limit n —> oo the mapping /(x) 
determines the state of the system, and particularly, 
the canonical ensemble averages of all physical quan¬ 
tities that we are interested in here. Consequently, 
Eqs. (10) — (13) along with Eqs. (16) and (17) provide the 
exact solution of the spin-1/2 Ising-Heisenberg model on 
a THE. 


C. Pairwise entanglement of the Heisenberg trimer 

In this subsection we turn to characterizing the entan¬ 
glement properties of the Heisenberg trimers of the THE. 
To quantify the amount of quantum correlations that the 
Heisenberg spins of each triangle share, we use below the 
concurrence [42] , as a computable entanglement measure 
of a bipartite mixed state. For this purpose we require 
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the reduced density matrix of a pair of Heisenberg spins 
(i.e., qubits), taking into account that due to the cyclic 
symmetry that each Heisenberg triangle possesses, the 
entanglement of any of its pairs is identical. We also note 
that as the neighboring Heisenberg triangles are coupled 
to one another by means of Ising-type (diagonal) inter¬ 
action, they are not entangled with each other. 

In order to express the reduced density matrix of a pair 
of triangular qubits (that we denote by P 12 ), in terms of 
the above defined magnetization toh smd the spin-spin 
correlation functions ^yy and we right it down 
in the Hilbert-Schmidt basis [48]: 

1 f ^ 

Pi2 = j (/(8)/-b2^ (ojCTi ®I + biI® ai) -I- 

\ i=i 

3 \ 

tyCTj 0 cTj- . (18) 

ib=i / 

Here / is a 2 x 2 identity matrix, and cTi’s are the standard 
Pauli matrices, related to the spin operators as = ici, 
= i(T 2 5^ = ^CTa. Following the well-known proce¬ 
dure we express the coefficients of the above decomposi¬ 
tion through the single-site magnetization and spin-spin 
correlation functions: 


0-3 = b3 = toh, til = 4 exx, 

t 22 = 4£yy, t 33 = (19) 

whereas all the other coefficients are equal to zero. As 
a result, the reduced density matrix of a pair of qubits 
reads: 
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f J+^Z 
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0 
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that is a special case of the so-called X state [49]. Its 
concurrence, as known, takes the following form: 


C(Pi 2 ) — 


2 max 2 \exx\ - 


( 21 ) 


+ £22 - mu 


1 


+ £22 + Toh , 0 


Note that this final expression for the pairwise entan¬ 
glement of the triangulated Husimi lattice is exact and 
takes into account the specific recursive structure of the 
system. 


III. BIFURCATION AND CHAOTIC 
BEHAVIOR OF THE MAGNETIZATION AND 
ENTANGLEMENT 


Having the exact solution of the spin-1/2 Ising- 
Heisenberg model on a THL we proceed to the discussion 
of its magnetic and entanglement properties. 


In what follows we are mainly interested in the case of 
isotropic antiferromagnetic Heisenberg interactions, i.e., 
= Jh < 0 , and the three-site interaction 
J 3 is also considered to be of antiferromagnetic charac¬ 
ter (J 3 < 0). Nevertheless, we choose the Heisenberg 
interactions to be dominant (Js/Jh < 1 ), which rein¬ 
forces our assumption of having an alternating sequence 
of Ising and Heisenberg variables (for weaker interactions 
the fluctuations in the transverse direction are expected 
to be reduced in the presence of a magnetic field di¬ 
rected along the z axis, which may give rise to spin- 
spin coupling of Ising, rather than of Heisenberg char¬ 
acter). Additionally, we assume the ferromagnetic Ising- 
to-Heisenberg coupling to be weak, by fixing the corre¬ 
sponding ratio to Jiu/Jn = — 0.01 (as already mentioned, 
the effective two-site Ising interaction J|® always remains 
ferromagnetic here). Furthermore, we set the magnetic 
fields acting upon the Ising and Heisenberg spins equal: 
Hu = Hi = H. Finally, we note that the initial [Eq. (1)], 
as well as the effective [Eq. ( 8 )] systems possess a particu¬ 
lar symmetry, namely, their properties remain unchanged 
under the replacement J 3 —>■ — J 3 and H —>• —H (below 
we consider the region of positive magnetic fields). 


A. Magnetization 

We start with the magnetic properties at relatively 
high temperatures, where the system possesses usual 
properties as shown in Eig. 2(a) — both for Ising- and 
Heisenberg-type spins the magnetization goes monoton- 
ically to saturation with the increase in the strength 
of the magnetic field. Nevertheless, as the tempera¬ 
ture is decreased we find the appearance of magneti¬ 
zation plateaus of the Heisenberg trimer at 1/6 of the 
saturation value, that corresponds to the ground states 

W = ^{\ m) + I nt) + I ;tt)) for i? > O and IV') = 
(I lit) + I iti) + I tii)) for i? < 0 (note that these 
states are highly entangled W states to which we return 
in the next subsection). 

The above picture, is, however, well known, and we 
do not go into more details in this respect. Instead, we 
consider the properties of the system at even lower tem¬ 
peratures, where, due to the specific lattice structure, it 
exhibits a doubling bifurcation — the mapping /(x) in¬ 
stead of a stable fixed point converges now to a stable 
two-periodic cycle. Consequently, the Ising mi magne¬ 
tization possesses two values, alternating from shell to 
shell, which can be interpreted as single-site magnetiza¬ 
tions of two emerging sublattices [Eig. 3 (a)]. Although 
this behavior has been studied quite in detail for a simple 
Ising model on a Husimi lattice [37, 39, 40], it results in 
a novel feature that we find here. Namely, the Heisen¬ 
berg Toh magnetization bifurcates as well, that is due to 
the consistency of the both subsystems. In other words, 
despite the ground state of the Heisenberg trimer being 
rigid, its environment (i.e., the Ising-type spins that con¬ 
stitute the Husimi lattice) brings about drastic changes 
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(a) HI\Jh 



FIG. 2. (Color online) (a) The single-site magnetization mn 
of the Heisenberg (red curve) and mi of the Ising (blue curve) 
sublattices versus the magnetic field 77 h = Hi = H . Here 
Jih/Jh = —0.01, Js/Jh = 0.5, T/IJhI = 1 and g = 3 (J^^ = 
Jh^ = Jh); (b) The same, but for T/|Jh| = 0.1. 


in its thermodynamic properties, resulting in a sublat¬ 
tice structure of the Heisenberg system, as shown in the 
inset of Fig. 3(a) (note also that this doubling bifurca¬ 
tion appears in the area of the formation of the above 
mentioned plateau at toh = 1/6). More interestingly, 
as we continue decreasing the temperature, we find fur¬ 
ther doubling bifurcations both for the Heisenberg and 
for the Ising magnetizations — new bubbles are formed 
as parts of the old ones, resulting in 2™ periodic phases 
(m = 2,3,...), with more complicated sublattice struc¬ 
ture. 

Finally, for ultimately low temperatures the system un¬ 
dergoes a transition to a chaotic phase, where the magne¬ 
tization pattern does not repeat itself anymore. We em¬ 
phasize again that the full period-doubling cascade and 
the chaotic behavior appears not only for the Ising mag¬ 
netization, but also for the Heisenberg one, which will 
consequently have a crucial impact on the entanglement 
behavior as well. Moreover, the chaos ensues in the vicin¬ 
ity of the incipient magnetization plateau at mu = 1/6, 
i.e., where the above mentioned highly entangled W state 
is the ground state of the Heisenberg trimer. 

Furthermore, the chaotic regime contains p-periodic 
windows (p = 3, 5,...) with p ■ 2™ periodic phases. Par¬ 
ticularly a wide three-periodic window is plainly distin¬ 
guishable in Fig. 3(b). Note also that the transition from 


and to chaos at both edges of the window takes place by 
means of a tangent bifurcation — behavior that is quite 
different from that of usual polynomial mappings, such 
as the logistic one. This peculiar feature of maps describ¬ 
ing spin-lattice models has been discussed in our previous 
works (note that the Q-state Potts model on a Bethe lat¬ 
tice exhibits similar behavior versus the temperature T 
and for noninteger values of Q < 2) [39, 40]. 




FIG. 3. (Color online) (a) The single-site magnetization mu 
of the Heisenberg (red curve) and mi of the Ising (blue curve) 
sublattices versus the magnetic field Hu = Hi = H. Here 
Jiu/Ju = —0.01, Js/Ju = 0.5, T/IJh] = 0.05 and g = 3 
{J^ = Ju = >/h). The inset shows the details in the area of 
a period doubling of the magnetization mu; (b) The same, but 
for T/\Ju\ = 0.02 (the inset shows the details of the chaotic 
behavior of the magnetization mu)- 


Meanwhile, to show explicitly the appearance of chaos, 
and to distinguish more rigorously between long-periodic 
phases and truly chaotic behavior, we consider the Lya¬ 
punov A(a;) exponent. The latter tells one whether an 
infinitesimal perturbation in initial conditions has an in¬ 
finitesimal effect (A(a;) < 0 — periodic behavior), or leads 
to a totally different trajectory (A(x) > 0), that would 
correspond to a chaotic regime (note that at a bifurca¬ 
tion point A(a;) = 0). The Lyapunov exponent is defined 
as: 


A(a;) 


lim — In 

m—¥co tTI 


dx 


( 22 ) 


where /^'"^(x) stands for the mth iteration of the map¬ 
ping f{x). 









































As Fig. 4 shows, the system indeed exhibits chaos, that 
corresponds to the areas where A(a;) > 0. Additionally, 
the above mentioned tangent bifurcations, corresponding 
to a transition from a chaotic phase to a periodic window, 
are also seen here [the periodic windows are the inclusions 
of negative Lyapunov exponent inside the chaotic phase; 
cf. Fig. 3(b)]. Finally, the model possesses a variety of 
superstable cycles with A(a;) = —oo [50]: as was shown 
in Ref. [39], at these points a specific change occurs in 
the symbolic dynamics of the mapping (16). 



FIG. 4. (Color online) The Lyapunov exponent A( 2 ;) of the 
mapping f{x) [Eq. (16)], versus the magnetic field Ha = Hi = 
H. Here Jih/Jh = —0.01, Js/Jh = 0.5, T/IJh] = 0.02 and 
q = 3 (Jh^ = Jh^ = Jh) [the values of the parameters are the 
same as in Fig. 3(b)]. 


B. Entanglement 

Before starting the analysis of the entanglement prop¬ 
erties of the Heisenberg trimer of the THL, we note that a 
system of three qubits has been already studied from dif¬ 
ferent perspectives and in various environments, includ¬ 
ing, e.g., external [51], or effective magnetic fields (that 
may describe a specific lattice structure [20]), incident 
[9] or cavity light field dressing [52], etc., and therefore 
is quite well understood by now. Nevertheless, the re¬ 
cursive structure, bringing about chaotic behavior that 
we present below, has not been considered yet. There¬ 
fore, we do not stop much in detail on the analysis of the 
ground-state features of the trimer of the THL, but pro¬ 
ceed, instead, to the studies of its periodic and chaotic 
entangled regimes. 

Having the exact formula for the concurrence C(pi 2 ) 
of a pair of qubits from the Heisenberg trimer, expressed 
in terms of the magnetization and spin-spin correlation 
functions [Eq. (21)], we plot it in Fig. 5, as a function 
of the magnetic field iLn = Hi = H. As Fig. 5(a) 
shows, the overall behavior of the concurrence at rela¬ 
tively high temperatures and weak Ising-to-Heisenberg 
coupling is quite analogous to that of a simple system of 
just three qubits in an external magnetic field: the figure 
compares the pairwise entanglement of a pure Heisen¬ 


berg trimer (dashed red curve), with that of our model 
(full blue curve). This similarity becomes evident from 
the fact that the coupling to the rest of the system, 
and particularly, to the Ising spins that constitute the 
Husimi lattice, happens by means of the weak Jih in¬ 
teraction. Meanwhile, as for the broad maximum at 
around C{pi 2 ) = 1/3, it appears due to the ground W 
state \iIj) = ^ (j tt/) + I tit) + I ttt)), that passes to 
a saturated non-entangled state | ttt) (the thermal ef¬ 
fects smooth here the abrupt jump from C(pi 2 ) = 1/3 
to C{pi 2 ) = 0). Note, however, that as Jm/Jn ratio 
is increased, the transition from the W to the saturated 
state is shifted to smaller values of H, since the Ising-to- 
Heisenberg coupling can be interpreted as an additional 
magnetic field of a strength JmPk {pk = ±3/2, ±1/2), 
acting upon the qubits [see also the Hamiltonian (1)]. 

Nevertheless, the features of our model at sufficiently 
low temperatures are strikingly different from the above 
picture. Namely, the entanglement of the Heisen¬ 
berg trimer bifurcates as the temperature is decreased, 
even when the Ising-to-Heisenberg coupling is weak [see 
Fig. 5(b), and its details in the area of period doubling 
in Fig. 5(c)]. This means that the concurrence C{pi 2 ) 
has now two values, that interchange one with another 
from shell to shell. In other words, we find here a rise 
of two sublattices with different values of pairwise entan¬ 
glement. Moreover, the breaking up into sublattices with 
respect to the entanglement is identical to that with re¬ 
spect to the magnetization mu: the upper branch of toh 
in Fig. 3(a) corresponds to the upper branch of C{pi 2 ) 
in Fig. 5(c). It is also worth noting that the bifurcation 
occurs in the weak magnetic field region, where the mag¬ 
netic entanglement is induced (the ground state at H = 0 
is not entangled, and quantum correlations build up here 
due to the inclusion of the magnetic field [53]). 

Meanwhile, at ultimately low temperatures, in the 
same manner as for the magnetizations mi and mu, the 
concurrence exhibits more and more phases with higher 
2™-periods (m = 2,3,...), and eventually, one reaches 
the regime of chaotic entanglement (Fig. 6). This means 
that in the entire system (deep inside the Husimi tree) 
the sequence of the values of the concurrence does not re¬ 
peat itself, and is, in other words, incommensurate. We 
also note that the above chaotic behavior appears solely 
due to the specific — recursive — structure of the THL, 
and due to the presence of three-site interactions in the 
initial (I), and therefore, in the effective Ising (8) model 
(as is known, the Ising model on a Husimi lattice with 
only two-site interactions does not exhibit chaos [37, 54]). 
Moreover, the region of the chaotic entanglement is di¬ 
rectly related to the strength of the three-site interac¬ 
tions — the stronger are these interactions, the wider is 
the chaotic region. 

On the other hand, we note that purely thermal entan¬ 
glement, i.e., entanglement that arises at H = 0, only due 
to the thermal mixing of the eigenstates of the system, 
exhibits neither period doubling, nor chaotic behavior. 
In other words, chaos ensues only in the presence of the 
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(a) til\JH\ 




HI\Jh\ 


FIG. 5. (Color online) The concurrence of a pair of qubits 
from a Heisenberg trimer versus the magnetic field. The 
full blue curves correspond to the concurrence C[pi 2 ) of our 
model, given by Eq. (21) (IFh = Hi = H , = Jh), 

while the dashed red curves stand for the concurrence of a sim¬ 
ple Heisenberg trimer with two-site isotropic interactions of a 
strength Jh, in an external magnetic held H (cf. [51]). Here 
Jih/Jh = —0.01, Js/Jh = 0.5, 5 = 3 and (a) T/IJh] = 0.09; 
(b) T/IJh] = 0.05; (c) Enlargement of the period doubling 
zone shown in (b). 


external magnetic field, where an entangled state appears 
in the ground-state structure. Furthermore, chaos does 
not induce entanglement by itself, but may enhance (or 
reduce) its amount from shell to shell. Specifically, as 
shown in Fig. 6(b), the chaotic entanglement is domi¬ 
nantly stronger than that of just three qubits in an ex¬ 
ternal magnetic field (compare the full blue and dashed 
red curves therein). 

Furthermore, an important observation is that the 
Lyapunov A(a;) exponent [Eq. (22)], that describes the 
system’s chaotic behavior, characterizes both the mag¬ 
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FIG. 6. (Color online) (a) The concurrence of a pair of qubits 
from a Heisenberg trimer versus the magnetic field. The 
full blue curve corresponds to the concurrence C{pi 2 ) of our 
model, given by Eq. (21) (iJn = Hi = H, = Jh), 

while the dashed red curve stands for the concurrence of a 
simple Heisenberg trimer with two-site isotropic interactions 
of a strength Jh, in an external magnetic field H. Here 
Jih/Jh = -0.01, Js/Jh = 0.5, T/\Jn\ = 0.02 and 5 = 3; 
(b) Enlargement of the chaotic zone shown in (a). 


netic and the entanglement quantities (e.g., the magneti¬ 
zation and the concurrence) of the present model. More 
precisely, as already mentioned in Sec. IIB, the prop¬ 
erties of the system in the thermodynamic limit are de¬ 
fined by the mapping (16), as particularly are the per-site 
magnetizations toh and mi, as well as the pairwise en¬ 
tanglement of the Heisenberg trimer, expressed in terms 
of < 7 (^ 12 ). Therefore the A(a;) exponent of the mapping 
(16), plotted in Fig. 4, confirms and quantifies the above 
discussed chaotic behavior of the spin-spin entanglement, 
shown in Fig. 6. Moreover, the inclusion of X{x) < 0 re¬ 
gions inside the chaotic regime, and A(a;) = —c» points 
in Fig. 4 reveal the existence of periodic windows and su¬ 
perstable cycles of the entanglement in the present model 
(similar to that of the magnetization, as discussed in the 
previous subsection). 

Finally, for more details, in Fig. 7 we show the transi¬ 
tion line between the phase without sublattice structure 
(uniform phase) and the 2-periodic one on the one hand, 
and the transition line between periodic and chaotic 
regimes on the other hand. More precisely, the upper line 
in Fig. 7 corresponds to the first bifurcation point, and 
up to the lower line (the shaded region) one has the full 
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2 ™ (m = 1,2,3,...) period doubling cascade with peri¬ 
odic entanglement behavior. Meanwhile below the lower 
line the system exhibits aperiodic chaotic entanglement, 
that, however, may be interrupted with p-periodic win¬ 
dows (p = 3, 5,...). 



FIG. 7. (Color online) The phase transition lines from the 
uniform phase to a two-periodic one (the upper curve) and 
from periodic to chaotic regimes (the lower curve). The full 
period-doubling cascade takes place in the shaded area be¬ 
tween these two curves. Below the lower curve the chaotic 
entanglement ensues, which may be interrupted with peri¬ 
odic windows. Here Jih/Jh = —0.01, Js/Jh = 0.5, and q — 3 
= jr = Jh). 


IV. CONCLUSION 

We proposed a multisite interaction spin-1/2 Ising- 
Heisenberg model on a triangulated Husimi lattice to 
study the effects of chaos and bifurcation on the system’s 
entanglement features. We used the generalized star- 
triangle transformation to map the initial model onto an 
effective Ising one on a simple Husimi lattice, which we 
then solved by means of the recursive method. An exact 
formula was obtained for the concurrence, to quantify the 
pairwise spin-spin entanglement of the embedded Heisen¬ 
berg triangles, expressed through the system’s single-site 
magnetization and spin-spin correlation functions. 

We have shown that at relatively high temperatures 
and weak Ising-to-Heisenberg coupling, the model ex¬ 
hibits quite usual magnetic behavior and its entangle¬ 
ment properties are analogous to that of just three in- 
tercoupled qubits. However, one finds drastic changes in 
the properties of the system as the temperature T is de¬ 
creased. Namely, period doubling of per-site Ising and 
Heisenberg magnetizations appears such that the corre¬ 
sponding values interchange one with another from shell 
to shell (this corresponds to a bifurcation point of the 
mapping that describes our model). We interpret this be¬ 
havior as a rise of a sublattice structure, that affects the 
system’s entanglement properties, too: the concurrence 
now exhibits period doubling as well, and each of the two 
values of the entanglement is related to a particular mag¬ 
netic sublattice. In other words, here every branch of the 


concurrence corresponds to a specific branch of the Ising 
and Heisenberg single-site magnetizations. 

On the other hand, with the further decrease in the 
temperature, we find a full-period doubling cascade, and 
eventually a transition to a chaotic phase where the se¬ 
quence of the entanglement (and magnetization) values 
is aperiodic and does not repeat itself. Importantly, this 
behavior appears due to the inclusion of three-site in¬ 
teractions — the system is always periodic when only 
pair interactions are taken into account. It is also in¬ 
teresting to note that chaos ensues in the region of in¬ 
cipient magnetic and concurrence plateaus, that corre¬ 
spond to the entangled W ground state. Additionally, 
chaotic behavior is absent if the system exhibits only 
thermal entanglement (in the absence of the magnetic 
field), i.e., if its ground state is separable, which, there¬ 
fore, along with the above argument, points to a con¬ 
nection of ground-state entanglement and chaos. More¬ 
over, although chaos does not induce entanglement by 
itself, it however may slightly enhance the latter (some¬ 
what similar results, namely entanglement enhancement 
for chaotic initial conditions, have been also indicated, 
e.g., in Ref. [31], although within distinct mechanisms 
for the chaos development). It is important to note, how¬ 
ever, that both for 2™-periodic (m = 1, 2,...) and for the 
chaotic regimes the ground state of the Heisenberg trimer 
is rigid, whereas the Ising sublattice, that is (weakly) cou¬ 
pled to the Heisenberg triangles, induces the underlying 
complex behavior. Consequently, at sufficiently low tem¬ 
peratures the entanglement (as well as the magnetiza¬ 
tion) of the Heisenberg spins exhibits a rich phase struc¬ 
ture, including the above mentioned 2™-periodic phases, 
chaos and p-periodic windows (p = 3,5,...), as well as 
various superstable cycles, that we have also confirmed 
by the studies of the Lyapunov A(a:) exponent. The 
above transitions between periodic phases occur through 
doubling bifurcations, whereas back and forth transitions 
from chaos to periodic windows happen through tangent 
bifurcations, that we plan to study in more detail in our 
future works. 

Finally, we emphasize that for the explicit studies of 
the chaotic behavior of the above spin model and specifi¬ 
cally, the chaotic pairwise entanglement of its Heisenberg 
trimer, we have mainly used the Lyapunov exponent ap¬ 
proach. The latter allows one to differentiate rigorously 
between periodic and chaotic regimes and to measure 
the strength of the emerging chaos. In this respect, our 
choice of the Lyapunov exponent characteristics has been 
particularly reinforced by the fact that the properties of 
the present spin model (and its state in the thermody¬ 
namic limit) are defined by a one-dimensional recursive 
mapping. As is known, in such a case the Lyapunov A(a;) 
exponent serves as a strong tool for the analysis of the 
inherent chaotic (as well as periodic) regimes [39, 40, 55], 
and moreover, for the model that we studied here, it char¬ 
acterizes the chaotic behavior not only of the system’s 
magnetic, but also of its entanglement features. Never¬ 
theless, it is worth noting that for a deeper understanding 
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of the relation between entanglement and chaos, apart 
from the Lyapunov exponent studies, one would necessi¬ 
tate here additional analysis, as, e.g., the estimation of 
the ergodicity of wave functions and level spacing statis¬ 
tics, quantum purity of the emerging mixed X states, etc. 
(see, e.g. Refs. [22, 23, 26]): an important issue that we 
will address in our future work. 
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